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We calculate the one-loop squared contributions to the next-to-next-to-leading order 0(a 2 af) 
radiative QCD corrections for the production of heavy quark pairs in the collisions of unpolarized 
on-shell photons. In particular, we present analytical results for the squared matrix elements that 
correspond to the product of the one-loop amplitudes. All results of the perturbative calculation 
are given in the dimensional regularization scheme. These results represent the Abelian part of 
the corresponding gluon-induced next-to-next-to-leading order cross section for heavy quark pair 
hadroproduction. 
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I. INTRODUCTION 

The increasing precision of present and forthcoming ex- 
periments in high energy physics requires a correspond- 
ing precision of the theoretical predictions. The next- 
to-leading order (NLO) predictions for hadronic heavy 
quark pair production suffer from inherent theoretical 
errors because of the well-known large uncertainty in 
choosing the renormalization and factorization scales. 
These errors are expected to be greatly reduced at next- 
to-next-to-leading order (NNLO) and therefore the need 
for a NNLO calculation of hadronic heavy quark pair 
production in QCD is by now clearly understood. 

A few years ago virtual two-loop and loop-by-loop cor- 
rections (we shall also refer to the "one-loop squared" 
contributions as the "loop-by-loop" contributions) were 
calculated by several groups in massless QCD (see e.g. 
[l| and references therein). The completion of a sim- 
ilar program for processes that involve massive quarks 
requires much more dedication and work since the inclu- 
sion of an additional mass scale dramatically complicates 
the whole situation. There are a number of publications 
where physicists work towards building up the necessary 
tools for calculating two-loop massive processes. For ex- 
ample, there are papers in which fully analytical forms of 
various master integrals are derived (see e.g. 0,0,1111). 
Recently Bernreuther et al. calculated the two-loop ver- 
tex corrections to heavy quark pair production from vec- 
tor and axial vector currents [g, Q. These results were 
utilized to determine a partial result on the forward- 
backward asymmetry in the process e + e" — > QQ involv- 
ing the two-loop contributions folded with the Born term 
and the loop-by- loop contributions Q. 

In general, there are four classes of contributions that 
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need to be calculated for the NNLO corrections to the 
hadronic production of heavy quark pairs. The first class 
involves the pure two-loop contribution which has to be 
folded with the leading order (LO) term. The second 
class of diagrams consist of the so-called loop-by-loop 
contribution arising from the product of one-loop virtual 
matrix elements which, for the special case of 77 colli- 
sions, form the subject of this paper. Further there are 
the one-loop gluon emission contributions that are folded 
with the one-gluon emission graphs. The one-loop gluon 
emission contributions also include the interesting class 
of the so-called pentagon graphs. Finally, there are the 
squared two gluon emission contributions that are purely 
of tree-type. 

In this paper, we concentrate on heavy quark pair pro- 
duction in photon-photon collisions which constitute the 
Abelian part of the gluon-induced hadroproduction of 
heavy quark pairs. From the technical point of view 
photon-induced heavy quark pair production is much 
simpler than the corresponding hadroproduction of heavy 
quark pairs. First of all there are no contributions from 
the subprocess qq — > QQ. Second there are no contri- 
butions from three-gluon coupling graphs which implies 
that there will be no collinear singularities and therefore 
the highest singularity in the one-loop amplitudes will be 
an infrared (IR) singularity proportional to (1/e). This 
in turn implies that the Laurent series expansion of the 
one-loop amplitudes can be truncated at 0(e) when cal- 
culating the loop-by-loop contributions. This is quite a 
bonus since it is the 0(e 2 ) terms in the Laurent series 
expansion that really complicate things || . Whereas the 
0(e ) contributions in the one-loop amplitudes involve a 
multitude of multiple polylogarithms of maximal weight 
and depth four 0, Q the 0(e) contributions needed in 
the present calculation involve at most trilog functions 
with their accompanying £(3) functions. 

It has been emphasized by many physicists that run- 
ning the ILC in the photon-photon mode is a very in- 
teresting option for the ILC (see e.g. [T3, EH)- The 
high energy photons can be generated by Compton back- 
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scattering of laser light on the high energy electron and 
positron bunches of the collider with practically no loss 
in energy and luminosity. Note that this reaction is also 
a clean channel for the investigation of various properties 
of heavy quarks (see e.g. [12] )• 

In this paper we report on a calculation of the NNLO 
squared one-loop matrix elements (loop-by-loop contri- 
bution) for the process 77 — > QQ. The calculation is 
carried out in the dimensional regularization scheme [l3j 
with space-time dimension n = 4 — 2e. The 0(e) expan- 
sion for the one-loop scalar master integrals that enter 
the calculation have been determined by us in For 
the divergent and finite terms the relevant amplitude ex- 
pressions were given in 1141 . The order e amplitudes have 
been written down in [151 ]. 

In a sequel to this paper we shall present results on the 
square of photoproduction amplitudes (75 — > QQ). Fur- 
ther we shall calculate the squares of the non-Abelian 
gluon-induced gg — > QQ amplitudes and the quark- 
induced qq — > QQ amplitudes which are needed for the 
loop-by-loop part of the NNLO description of heavy fla- 
vor production. 

In our presentation we shall make use of our nota- 
tion for the coefficient functions of the relevant scalar 
one-loop integrals calculated up to 0(e 2 ) in For the 
77 — * QQ case one needs one scalar one-point function 
A, the four scalar two-point functions Si, S2, S3 and S4, 
the three scalar three-point functions C2 , C5 and Cq , and 
one scalar four-point function D\ . As was mentioned be- 
fore, the 77 — * QQ amplitudes are only IR-singular due 
to the absence of a three-gluon coupling contribution. 
And, in fact, one finds that the above set of multi-point 
amplitudes have at most a 1/e singularity [5j . For exam- 
ple, employing the notation of [5j, we define successive 
(i) 

coefficient functions D[ for the Laurent series expan- 
sion of the complex scalar four-point function D±. One 
has 

Si = zC £ (m 2 ){^D[- 1 '> +eS« + 0{e 2 )), (1.1) 

where C £ (m 2 ) is defined by 

We use this notation for both the real and the imaginary 
part of Si, i.e. for KeDi and ImDi. As also mentioned 
before we can truncate the Laurent series at 0(e) since 
this is sufficient for the 77 — > QQ case treated in this 
paper. Expansions similar to Eq. (jl.ip hold for the scalar 
one-point function A, the scalar two-point functions Bi 
and the scalar three-point functions Ci. 

The paper is organized as follows. Section HT1 contains 
an outline of our general approach as well as a discussion 
of the singularity structure of the NNLO squared matrix 
element. In Section Hill we discuss the structure of the 
finite part of our result. Our results are summarized in 
Section IIVI Finally, in an Appendix we present results 



for the various coefficient functions that appear in the 
main text. 



II. NOTATION AND THE SINGULARITY 
STRUCTURE OF THE SQUARED AMPLITUDES 

The one-loop Feynman diagrams relevant for heavy 
flavor production by two on-shell photons 

1(P1)+7(P2)^Q(P3) + Q(P4), (2.1) 

are depicted in Fig. [TJ 

The directions of the momenta correspond to the phys- 
ical configuration, e.g. p\ and P2 are ingoing whereas p 3 
and P4 are outgoing. With to the heavy quark mass we 
define: 

s=( Pl +p 2 ) 2 , t = T -to 2 = (pi -p 3 ) 2 -to 2 , 

u = U - m 2 = (p 2 -p 3 ) 2 - m 2 . (2.2) 

so that the energy-momentum conservation reads s + 1 + 
u = 0. 

When squaring the amplitudes one sums over the spins 
of the final state heavy quarks. We have decided to also 
average over the spins in the initial state when presenting 
our results. In the spin-averaging over the initial photons 
we divide by a factor of four. One could as well divide by 
(n— 2)(n— 2) where (n — 2) — 2(1— e) are the spin degrees 
of freedom of a on-shell massless photon in DREG. This 
is a matter of convention. If one wants to convert to the 
convention where the spin-averaging factor is (n— 2)(n — 
2) one has to multiply our result by the overall factor 

(1-e)" 2 = l + 2e + 3e 2 + C(e 3 ). (2.3) 

Also, our results have to be multiplied by the overall 
factor 

C = (g 2 e 2 Q g 2 C £ ( m 2 )) 2 , (2.4) 

where g and g s are the renormalized electromagnetic and 
strong coupling constants, respectively, and eq is the 
fractional charge of the heavy quark. C e (m 2 ) is defined 
in fL2]) . 

When squaring amplitudes, we contract polarization 
vectors in the Feynman gauge. Note that we have used 
the on-shell conditions»i • e\ = and pi ■ £2 — in our 
amplitude calculation [Ij,[lB] according to the framework 
set up in [l6j]. This has advantages in the non-Abelian 
case since one can omit ghost contributions when squar- 
ing the amplitudes. Using the above on-shell conditions 
already on the amplitude level means that one takes full 
advantage of gauge invariance when squaring the ampli- 
tudes. 

As shown e.g. in [13, [3 the self-energy and vertex di- 
agrams contain ultraviolet (UV) and IR poles after mass 
renormalization. It is well known that the renormaliza- 
tion of the wave function and QED-type vertex renormal- 
ization effectively imply multiplication of all the external 
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FIG. 1: The f-channel one-loop graphs contributing to the photon fusion amplitude. Wavy lines represent the photons, curly 
lines represent the gluons and the thick solid lines correspond to the heavy quarks. 



massive quark self-energy graphs by one-half. We used 
this fact as a partial check of our calculation. Halving the 
external self-energies (in our case these are matrix ele- 
ments for the graphs Fig.[ljd2) and Fig.[TJd3)) we have 
explicitly verified cancellation of UV divergences for all 
the relevant sets of squared amplitudes. 

In order to fix our normalization we write down the 
differential cross section for 77 — > QQ in terms of our 
amplitudes squared \M\ 2 . One has 



da 



77->QQ 

dtdu 



1 rf(PS); 

2 s dtdu 



(2.5) 



where the n-dimensional two-body phase space is given 
by 



d(PS) 2 



(4tt) e 



8tts r(l-e) 



tu — smS 



2\ -e 



5(s + t + u)dtdu. 

(2.6) 



The final spin summed and initial spin averaged 
squared matrix element can be written down as a sum 
of three terms: 



|M' 2 



Loop x Loop 



CO 



^ V (-2) + l V i-i) + NcV (^y 



(2.7) 

where Cp = 4/3 and C has been defined in (|2.4[) . Note 
that we used 1/4 for the spin-averaging factor as ex- 
plained before (|2.3p . The first two terms contain IR poles 
of the second and the first order, respectively. The third 
term represents the finite part. All three terms are origi- 
nally bilinear forms in the coefficient functions that define 



the Laurent series expansion of the scalar integrals (|1.1[) . 
Some of these coefficient functions are zero and some of 
them are just numbers. In the latter case we have substi- 
tuted these numbers for the coefficient functions in the 
three terms above. This has been done for the coefficient 
functions A«, B { ~ 1] , B { 2 ^\ B ( 3 j) and b[ j) . 

The most singular first term in (|2.7p is proportional to 
the leading order term of the process 77 — > QQ, e.g. 



iyHO = 1 
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(ge Q )' 



\M- 
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where the Taylor series expansion of the square of the 
LO amplitude |Mlo| 2 can be calculated to be 
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}. (2.9) 
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Note that the original amplitude expression contains the 
coefficient functions £)j ^ and d[ 
been substituted for by the (t «-> u 



(-1)1 



which have 
symmetric coeffi- 
cient function Cr 1 ' using the relations d[ ^ = Cg ^/t 



and D 



(-1)1 



= Cg ^ /u [|[. It is only after these sub- 
stitutions that the factorization of |Mlo| 2 in (|2.8p be- 
comes apparent. We remind the reader that the coef- 
ficient functions are complex functions. This has to be 
taken into account when calculating the three contribu- 
tions in Eq. (j!T7]) . 

For the second term in (|2.7p we obtain: 



J 



,{((2, 



s)C t 



(-1) 



- 1 



(2.10) 



-Cfh c o + |cf /i c ° + \cf ] h c o + D^h D o + (i «- u)] *} 



We have chosen to label the coefficient functions hi ac- order to obtain (|2.10p in a factorized form, 
cording to the coefficients of the Laurent series expansion 
which they multiply. Note we have again reexpressed the 

coefficient functions D^T 1 ' and D\ l '\ via Ci -1 " 1 in The relevant coefficients read: 



h = 2m 4 (16m 4 s + 4m 2 si + 28m 2 su - 8m 2 tu + 2stu - 19tu 2 - 30m 2 tu 2 /T)/(Utu 3 ) 

-(56m 8 u + 29m 6 st - l7m 4 t 2 U - 105m 4 tuT + 7m 2 t 3 T + m 2 t 2 uT - 2m 2 t 2 uU - 3st 3 u)/(s{3 2 tuTU), 
h B o = (4m 4 s 2 - 16m 4 sT + 8m 2 tTu - 3st 2 {2m 2 ~u)-t 4 (2m 2 - u)/T)/(t 3 u), 

h B o = -AD /(tup 2 ), (2.11) 

h c o = 2(6m 4 si - 2m 4 sw + 8m 4 tu - 2m 2 st 2 - 4m 2 t 3 - 2st 3 + t 2 u 2 )/(t 2 u), 

h c o = -s(8m 4 - 3s 2 + 2tu)/(tu), 

h c o = -(8m 4 s + 2m 2 s 2 + 8m 2 tu - s 3 - 2stu)/(tu), 

h D o = (16m 6 s + 24m 4 st - 8m 4 u 2 + 2to 2 s<m + 2s 2 t 2 - stu 2 )/{tu), 



where (3 = yl — 4m 2 / s. The terms multiplied by 1/2 
in (|2.10[) are (t <-> u) symmetric which, when adding the 
(t <-> u) term, add up to the full contribution. 

We emphasize that the real part of the second factor in 
square brackets in (I2.10| is nothing but the finite part of 
the NLO contribution (up to an overall trivial factor, see 
below) calculated in [17[ by summing the two finite terms 
(i.e. not multiplying the IR pole 1/e) in Eqs. (16) and 
(22) of [13] ■ We mention that one can also obtain the 
imaginary part of the second factor in square brackets 
from [ijj by using the procedure described in the last 
section of |14[- When comparing with the results in [l7| 
one has to of course substitute the explicit expressions 
for the coefficient functions given in []|. 

In the following we want to exhibit a remarkable struc- 
ture of the NLO and NNLO single pole contributions. For 
the NLO case we have 

\M^o\ 2 = - £ 9 2 s C £ (m 2 )C F (2.12) 

x Re {(2m 2 - s) - l} |A^lo| 2 ; 

Let us now write Eq. (|2.10p in the form 

= jN c (ReA ReB + ImA ImB), (2.13) 

where A and B stand for the round and square brackets, 
respectively, i.e. 

A = (2m 2 - s)C { f 1} - 1, (2.14) 

B = h + Bf ] h B o + \B { ° ] h B o + dphci 



+~CfV° + \c^h c « + D^h D o + (t~ u). 
One observes the remarkable pattern 

|MP°V = CC 2 -V^ = - £ g 2 C £ (m 2 )C F (2.15) 

x {ReA \M^l% c \ 2 + A(ge Q ) 4 g 2 C s {m 2 )C F N c lmA ImB} . 

We also mention that when calculating separately the 
t—,u— channel contributions and their interference the 
second power of the denominator D = m 2 s — tu would 
still be present in the results for the 1/e contributions. 
Only after adding up all contributions the two remain- 
ing powers of this denominator for the 1 / e contributions 
cancel. It is in this way that the factorized structure of 
the amplitude is obtained. In the previous NLO calcula- 
tion by one of us [13 , the cancellation of all the powers 
of the denominator D as well as the factorized structure 
was obtained separately for the t ® t or t ® u contribu- 
tions. Of course, when adding up all the channels, we 
obtain the same result as in [l?} for the complete gauge 
invariant set of diagrams. 

Note that the factorized form of Eq. (|2.15p holds only 
when one retains the full e dependence in the Born term 
as given in (|2.9p . 

We would like to mention that the nice and simple fac- 
torized form of (|2.15[) could have been anticipated from 
the fact that in the one-loop matrix elements the sin- 
gle pole always multiplies the LO term. Consider e.g. 
Eq. (7) of the recent publication [l8| where the predic- 
tions for an overall singular structure of massless one-loop 
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amplitudes are extended to the case of massive partons. 
In our case with only two colored massive particles, the 
second term in Eq. (9) of fl8j , i.e. containing the infrared 
factor which controls color correlations, completely fac- 
torizes. The second single pole structure comes from the 
function of Eq. (18) of Together, they give a full IR 
structure that is proportional to 

2C F 1 2m g ~ S (bxx + wr) - 1 j . (2.16) 

Considering the fact that the 1/e pole in our integral 
function Cq is multiplied by (In x + iir) / s/3, we indeed re- 



produce in our Eq. (|2.15[) the infrared structure derived 
(e.g. by formal squaring of Eq. (7)) from general con- 
siderations in 18]. In the next section we present the 
finite parts of our squared amplitudes. Quite naturally, 
these results cannot be derived from general principles 
and must be obtained from an explicit calculation. 



III. STRUCTURE OF THE FINITE PART 

In this section we present the finite part of our result. 
It can be written as 



y(o) = 1/ + Re {s(°) [/ B? + B^f B o B o + \B^f B o B o u + B ( pf B o B o + ci 0) f B o c o + c!?Jf B » c o u (3.1) 

1 f I 1 r>(°) f 
2JB" + 2^2 JB°B° 



+ G 5 JB°C° + ^6 JB°C° + V l JB°D° + U lu J B®D® u 



B?f B] 



B. 



(0) 



+ C 2 JB°C° + 2°5 JB°C° + 2°6 JB°C° + U l J B°,D° 



1 R (l) f 

1 U 2 JE 

.(1) 



2^2 3B\ 



a 



(0) 



fc° + C 2 0) fc°C° + h C 2ufc°C° 



(0). 



5^5 



(0) 



fc° + Cf ] f c o c a + C ( 6 0) f c o c o 



+2D { °>f c o D o 



+ ir (1) f , + ir (0) 

+ 2°5 /CJ + 
>(0) 



2^5 t- 2^6 



/eg + ^6 fc°C° + 2D ^f c o D 



1^ 

-.(0). 



(-i) 



D7 L B> 



12 12 15 15 lb 



+^6 /ort ■ 



^ 1 +-D 1 0) /n-nO+M 1) /n- lD J 



D 



(-1) 



D 



(0) 



R (0) f , R (l) f 



+ °2 lD7 x C° 



I 



where the subscript "u" is an operational definition pre- 
scribing a (t u) interchange in the argument of that 



function, i.e. B 



(0) 



B 



(0)1 



etc... Note that all three 



functions C, 



(-D n(-D 



L ^ and d[ u ^ appear when calculat- 



ing (|3.ip . as well as some spurious poles 1/e multiply- 
ing them. However, when substituting Cg in terms of 
D ( f 1] and I^ -1) | t „ u these spurious poles cancel out. 

Again the terms multiplied by 1/2 in (|3.ip are (t <-> it) 
symmetric which, when adding the (t <-> u) term indi- 
cated at the end of (|3.ip . add up to the full contribution. 

Several coefficient functions in (|3.1[) are trivially re- 
lated to the corresponding ones in (|2. 10p : 



h 



b" ■ 



fc 



-8/1/70 , 



(3.2) 



/ci = -8h c o, f c i = -8h c o, f D i = -8h D o. 



(cf. Eq. (jl.ip ) on the one hand, and powers of e resulting 
from doing the spin algebra in DREG on the other hand. 
For example, one has a contribution from D^^'B^ as 
well as a contribution from d[ B± . Terms of the type 

D^'B^, where the superscripts corresponding to e- 
powers do not compensate, would be absent in the reg- 
ularization schemes where traces are effectively taken in 
four dimensions, i.e. in the so-called four-dimensional 
schemes or in Dimensional Reduction. 

Finally we note that obtaining our results as well as 
casting them into the above compact forms (also for the 
coefficient functions presented in the Appendix) was done 
with the help of the REDUCE Computer Algebra System 



IV. 



CONCLUSIONS 



The results for the other coefficient functions appearing 
in (|3.ip can be found in the Appendix. 

In the finite contribution Eq. (|3.ip one can see the 
interplay of the product of powers of e resulting from 
the Laurent series expansion of the scalar integrals 



We have presented 0{a 2 a 2 s ) NNLO analytical results 
for the loop-by-loop contributions for heavy quark pair 
production in 77-collisions. The present paper deals 
with unpolarized photons in the initial state. Using the 
backscattering technique it is not difficult to obtain po- 
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larized photon beams of high intensity at the 77 option 
of the ILC by colliding the low energy laser light with 
polarized electron and positron beams. Since our ampli- 
tude results [IH contain the full spin information of the 
process an extension of the present paper to the case of 
polarized 77 production of heavy quark pairs would not 
be very difficult. 

The present results form an Abelian subset of the 
non-Abelian gluon-induced loop-by-loop contributions 
to heavy quark pair production. This calculation con- 
stitutes a first step in obtaining the analytical results for 
the exact NNLO corrections to the heavy quark produc- 
tion processes in QCD. Analytical results in electronic 
format for all the terms in Eq. (|2.7p are readily available 
[20l |. The next step would be to provide similar results 
for the two remaining channels of heavy flavor produc- 
tion that were discussed in the Introduction. We reserve 
this task for future work. 
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APPENDIX 

In this appendix we present the coefficient functions 
for the finite part of the loop-by-loop contributions ap- 
pearing in Eq. (|3.1[) . 

In order to keep our expressions as compact as possible 
we introduce the notation: 



zt — 2m 2 +t, z u = 2m 2 + u, D = m s — tu. 

(A.l) 

The coefficients of the finite part read: 



/ = -2{D0 2 TU2m 2 s{l2d>m w {t 2 + u 2 ) 2 - 208m 8 s(t 2 + u 2 ) 2 + 8m 6 (10i 6 + 39t 5 u + 34t 4 u 2 + lUt 3 u 3 + 34tV 
+39tu 5 + 10m 6 ) - 8m 4 stu(lM 4 + U 3 u + 66tV + Atu 3 + 13m 4 ) + m 2 t 2 u 2 {36t 4 + 7ht 3 u 
+190t 2 u 2 + 75tu 3 + 36m 4 ) - st 3 u 3 (7s 2 + 8tu)) 
+DTU2m 2 t 2 u 2 (12m 2 (t 5 + m 5 ) + tu{U 4 + 7t 3 u + 42< 2 u 2 + 7tu 3 + 4u 4 )) 

+s[3 2 TU'i8Am 8 t 4 u 4 + TUt 4 u 4 (2(t 4 T + u 4 U) - 4tu(t 2 T + u 2 U) - 6m 2 tu(t 2 + u 2 ) - 12iV(T + U) 
-s 5 + 12siV) - t 4 u 4 {t - u)(t 4 U(M - u) + m 4 T(< - 3u)))/(D{3 2 T 2 U 2 st 4 u 4 ), 

f B o = 4:(-T 2 UD(512m s T(t 2 + u 2 ) + 448m 8 (< 3 + u 3 ) - 704m 8 stu + 464m 6 i 4 + 1776m 6 t 3 u + 1216m 6 iV 

+688m 6 tM 3 + 80mV + 184m 4 i 4 M + 768m 4 i 3 u 2 + 288m 4 i 2 u 3 + 72m 4 tM 4 + 30m 2 iV - 23Tt 3 u 3 

+2m 2 t 2 u 4 - 22t 3 u 4 ) 
+T 2 Um 4 tu{7AQm 2 t 4 + lUm 2 t 2 u 2 + 12m 2 u 4 - 31i 3 u 2 + 33t 2 M 3 ) 
-T 2 p 2 st 2 (32m 8 t 3 + 32m e u 3 U + 12m 4 t 3 u 2 - 8m 4 su 4 - 16t 3 M 4 ) 
+T 2 t 4 u 3 {24 r lm 4 t - At 3 + 3u 3 ) - TDm 4 t h {A8m 2 t + 174m 2 u + 190tu) 

-T2t 5 u 2 (18m 4 t 2 - m 4 tu + 2m 2 t 2 u + m 2 u 3 - 2t 3 u) + Dt 7 u 2 (3m 2 + 4u) + m 2 t 7 u 4 )/(Df3 2 T 2 Ust 4 u 2 ), 

/b°b° = -2(DT(64m 10 + 8m 8 (29t + 5m) + 304m 6 i 2 + 96m 6 tu + 136m 4 i 3 + 68m 4 i 2 M - 12t 5 - 9t 4 u) 
+Tt 2 (2m 6 u 2 + 22m 4 t 2 T + m 2 t 2 (t 2 + 30tT + 5tu - 2m 2 ) + 5< 5 ) + t 8 )/(DT 2 t 4 ), 

/bJsJ = 4(-D(64to 12 - 168m 10 s + 8m 8 (19s 2 + 22tu) + 40m 6 (i 3 + u 3 ) - 398m 6 stu + 98m 4 tu{t 2 + u 2 ) 
+286m 4 tV - 49m 2 siV + 2t 3 u 3 ) + 6m 4 tu(2m 4 (t 2 + m 2 ) - t 2 u 2 ))/(DTUt 2 u 2 ), 

f B o B o = -8z t (2D(m 2 z t + t 2 ) - M 2 (4m 2 T + t 2 ))/(D(3 2 Tt 2 ), (A.2) 

f B o c o = 8(D 2 T2(6m 6 s ~ 2m 4 t(5t + 2u) - t 2 z t {At + u)) - DTt 2 (2z t (8m 6 - 7m 2 t 2 - 4i 3 ) - i 2 M 2 ) 
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+Dm 2 t 5 u ~ 3t 6 zf)/(D 2 Tt 3 ), 

f B °c° u = 8(D 2 Tm 2 (l2m 4 s - 4m 2 (3t 2 + 7tu + 3u 2 ) - u(l5t 2 +tu + 8m 2 )) 

-DTm 2 u(16m 4 uz u + m 2 t(YAt 2 + 31?i 2 ) + Wt 3 u + 10tu 3 ) + DU2t 2 u(m 2 t 2 + 3raV + tu 2 ) 
+Dm 2 t 2 u 2 (7m 2 u + t 2 ) - t 4 u 2 (2z t {<6m 2 U + u 2 ) - tuz u )) / [D 2 Tt 2 u), 

f B o c o = A{D 2 T{8m 4 {ht + it) - Am 2 s(t + 3u) - Atist + 2m 2 )) + D 2 t 2 (Ut 2 + lOtu + I3u 2 ) 

1 5 

-Zm 2 (8m 4 (lli + 19m) + 48m 2 t 2 + lit 3 ) - D(47m 2 t 4 u + Ylm 2 t z u 2 + 4m 2 t 2 u 3 - t 5 u) 
+T(32m 6 t 2 u 2 + 32m 4 i 2 u 3 + 12m 2 t 4 (s 2 - 3tu) + 3t 7 ) - 32m 12 s 2 + 3m 2 t 4 u 2 (4ra 2 + u))/{D 2 Tt 2 ), 

f B o c o = 4/3 2 s(L>(8m 8 s + 8mhu + 8m 4 t 2 (3t + Au) + m 2 t 2 (3lt 2 + 28tu - 9u 2 ) - t 3 u 2 ) 
-Am 8 s{t 2 - u 2 ) + 4m 2 i 3 (3si 2 - 2t 2 T - u 3 ) + st 5 (3t + 5u))/(D 2 Tt 2 ), 

fB°D° = 4(L» 2 T16m 6 (llt + 3m) + DT(64m 10 s - 4m 6 (llt 3 + I0t 2 u + Yltu 2 + 2u 3 ) 

+2m 4 t(51t 3 - At 2 u - 33tu 2 - 6u 3 ) + 2m 2 t 3 (31t 2 + 32tu + 7u 2 ) + t 3 (l8t 3 + \ht 2 u + 10tw 2 + 3u 3 )) 
-Dt 5 (7t 2 + su) - 3t 6 (16m 4 T + m 2 t(6t - u) + t 2 T))/(D 2 Tt 2 ), 

f B o D o u = 4(L> 2 T4m 2 (16m 6 + 4m 4 (7i + 3u) + m 2 (16t 2 + 27tu + u 2 ) + t 3 ) 

-L>T(64toV + 2m 4 u(23t 3 - 6t 2 u + 15tu 2 - 4u 3 ) + m 2 tu(6t 3 + lbt 2 u + tu 2 - 2u 3 ) + st 2 u 2 (t + 2u)) 
-Dm 4 t{At 4 + 3t 3 u - 26*V - \3tu 3 + 4u 4 ) + 3mH 3 u(T(t 2 + 15u 2 ) + u(Ast - tu + u 2 )))/(D 2 Tt 2 ), 

f B o = 8{8m 2 (2m 2 s 2 - 2m 2 tu - stu)/{t 2 u 2 ) + (24m 6 sf3 2 /D - 28m 2 sf3 2 - 11D + 3tu) / (TU)) / 1 , 

/ s o B o = 4(4m 2 s - 3t 2 + 2tu - 3u 2 )/(D(3 2 ), 

f B o c o = -8(L»2(10m 4 t + 2m 4 u + 3m 2 t 2 + m 2 tu + 2t 2 u) + 2 st 3 (lOm 2 + 3t) - tu(4m 2 uz t + t 2 (t - u) j)/(D 2 f3 2 ), 

f B o c o = As{DA{3m 2 s - t 2 - u 2 ) - 6m 2 s(t 2 + u 2 ) + 3{t 2 + u 2 ) 2 )/D 2 , 

f B o c o = 4s(L>4(3m 2 s - (t - uf) + 3t 3 (4m 2 +t) + 3w 3 (4m 2 + u) + 6t 2 u 2 )/D 2 , 

Ib°d° = 4(L> 2 8m 2 z„ + D4m 2 t 2 (5sf3 2 + 2t) - (3 2 6m 2 st 2 (3t 2 + u 2 ) + I0m 2 s 2 t 3 - st 3 (3s 2 - 2tu))/(D 2 f3 2 ), 

f c o = -8{D 2 {8m w u{A2t{t 2 + u 2 ) + t 2 u + 15m 3 ) + 2m 8 w(149t 4 + 595t 3 u + 95t 2 u 2 + 109tu 3 + I2u 4 ) 

+2m 6 t(4t 5 + 35t 4 u + 159t 3 u 2 + 328t 2 u 3 + 81tu 4 + I7u 5 ) 
+2m 4 t 2 u(25t 4 + 322t 3 u + \32t 2 u 2 + 7htu 3 + lOu 4 ) 

+2m 2 t 3 u(7t 4 + I3t 3 u + 9t 2 u 2 + 46tu 3 + 8u 4 ) + t 4 u 2 (10t 3 + 20t 2 u + \7tu 2 + 7u 3 )) 
+Df3 2 2m 2 st 2 (40m 8 u 3 + I00m 6 t 2 u 2 + 8m 4 t 5 + 2Alm 4 t 3 u 2 + 67m 4 tu 4 + 35™ 2 t 3 u 3 + 29t b u 2 + 62t 3 u 4 ) 
-Dm 2 (32m 12 u(t 2 + u 2 )(t + 3u) - 128m 10 t 5 - t 5 u 2 (49t 3 - 33t 2 u + 118tu 2 + 66u 3 )) 
+p 2 8m 2 st(5m 12 (t 4 - u 4 ) + 2m w t 3 (t 2 - llu 2 ) - 6m s tu 5 - t 3 w 4 (16m 6 + 43m 4 t + llm 2 t 2 + 9t 3 )) 
-48m 2 t 6 u 4 (t - u) 2 )/(D 2 f3 2 TUst 3 u 2 ), 

fc°C° = 4(L> 2 (4m 8 .s 2 + t 2 (l6m 6 t + 50mH 2 - 2m 4 u 2 - 28m 2 t 3 - 8t 4 - 2t 2 u 2 )) 
+Df3 2 st 2 (34m 6 t 2 - 2m 6 u 2 + 20m 4 t 3 + 2m 4 stu - 8m 2 t 4 - 4t 5 ) 
+D4t 5 (m 2 t 2 + 3m 2 u 2 + 2tu 2 ) - p 4 3s 2 t 8 )/(D 3 t 2 ), 



8 



/ c o c o = 8(D 3 A(m 6 s - m 4 (s 2 + 7tu) + tu(s 2 - tu)) + D 2 t 2 u 2 (76m 4 + (t - u) 2 ) 
-D2m 2 t 2 u 2 (6Am 6 + Am 2 tu - bstu) - 3m 4 t 2 u 2 {t - u) 4 )/(D 3 tu), 

f c o c o = A{D 3 8m 4 {7t + it) - D 2 Am 2 t{35m 2 t 2 + m 2 u 2 + 9t 3 - 2At 2 u - Atu 2 - u 3 ) 

-DT8m 2 t 2 {A0m 4 s + llm 2 t 2 - 24raV + tu 2 ) 
-D2m 2 t 2 {Am 4 {m 2 + 8u 2 ) + t 3 (A5t + 59u) - u 2 (5t 2 - tu - Au 2 )) 

+2m 2 t 3 (Ut 5 + I0t 4 u + A2t 3 u 2 + 15iV + I8tu 4 + u 5 ) - st 4 {3s 4 - 2t 3 u - 2su 3 ))/{D 3 t) 1 

f c o c o = A0 2 s(D 2 A{m 4 {9t + u) + m 2 (4i 2 - u 2 )) + DTA8t 2 (m 2 u + t 2 ) 

-D2{m 2 t(bt(t 2 + u 2 ) + 2su 2 ) - t 2 (Ust 2 + t 3 - lOstu + su 2 )) - 3t 3 {t - u) 4 )/D 3 , 

f c o D o = 4(£> 2 (16m 8 (3i 2 - 2tu - u 2 ) + 4m 6 (5i 3 - I2t 2 u - \\tu 2 - 2u 3 ) + 8m 4 t{8t 3 - 8t 2 u - 5tu 2 - u 3 ) 

-Am 2 t 2 (29t 3 + \2t 2 u + 7tu 2 + u 3 ) - 2t 3 (2t 3 + VM 2 u + htu 2 + u 3 )) 
-D2t 4 {128m 8 + 32m e t - 2AmH 2 + 3m 2 t 2 u - \7t 3 T) - M 7 {16m 2 tz u + (t - u)(s 2 - 8tu)))/{D 3 t), 

fc°D° u = 4(D 2 (16m 8 (3i 2 - 2tu - u 2 ) + 4m 6 (10t 3 - I3t 2 u - u 3 ) - Am 4 t(2t 3 + 3t 2 u - 16tu 2 + 3u 3 ) 
+Am 2 (st 2 (2s 2 - tu) - tu 4 )) 
-Dl&m 4 t 2 u{8m 2 uz t + (t - u)(t 2 + u 2 )) - 6m 4 t 2 u 2 (5t 4 + lOiV + u 4 ) - 3s 3 t 3 u 2 (2m 2 s - tu))/(D 3 t), 

f c o = A(D 2 (6Am la (t 2 + u 2 ) - m 6 (40s 4 + 96s 2 tu + \AAt 2 u 2 ) + m 4 (2As 5 + AAtu(t 3 + u 3 ) + 548s*V) 
-Am 2 tu(3s 4 + 55tu(t 2 + u 2 ) + 25* V) + st 2 u 2 (l7s 2 - Utu)) 
-Dm 2 ((A8m s - 25*V)(* - u) 4 - 32m 2 *V(2™ 2 - s)(8m 4 - 7tu) + t 3 u 3 (21s 2 - 116tu)) 
-2Am 2 t 4 u 4 (8m 4 + s 2 - 6tu))/(D 2 TUt 2 u 2 ), 

f c o c o = ~(D 3 896m 4 - D 2 2(8m 4 (27s 2 - 80tu) + 3Am 2 s 3 + 2Am 2 stu - A9(t 4 + u 4 ) - IQsHu + 2U 2 u 2 ) 

5 5 

+DA(128m 4 t 2 u 2 - 7(t 5 - u 5 )(t - u) + 12tu(t 4 + u 4 - t 2 u 2 ) - 22t 2 u 2 (t 2 + u 2 )) + 3(t 2 + u 2 ) 2 (t - u) 4 )/D 3 , 
f c o c o = 2p 2 s 2 (D2(A0m 6 s - 2m 4 (3s 2 - Atu) - 3m 2 s(7s 2 - 16tu) + lis 4 - £w(49s 2 - 50tu)) - 3(t 2 + u 2 )(t - u) 4 )/D 3 , 
f c o D o = 2(L> 3 64m 6 - £> 2 (32m 6 (18t 2 - Mu + u 2 ) - 24m 4 (13t 3 - 7t 2 u + btu 2 + u 3 ) 

5 1 

-Am 2 (A7t 4 + 77t 3 u + 3At 2 u 2 + lltu 3 + 3u 4 ) + 2t(2t 4 - t 3 u - 3t 2 u 2 - 21tu 3 - u 4 )) 
+D2t 2 (2b6m 6 t 2 + 288m 4 t 2 u - 8m 2 (9t 4 - tu 3 + u 4 ) - t 3 {\7t 2 + 59u 2 )) - 3t 2 {t - uf{uz t - t 2 ))/D 3 , 

f c o = A{-D 2 (8m e (3s 4 - 8s 2 tu + 22t 2 u 2 ) + Am 4 stu(22s 2 - U9tu) + Am 2 t 3 u 3 - st 2 u 2 (As 2 + 33tw)) 

+D(6Am 12 s(t 2 + u 2 ) - 16m 10 (3t 4 + 2t 2 u 2 + 3m 4 ) + 336m 6 tV(tT + uU) - 298m 4 st 2 u 2 (t 2 + u 2 ) 

+Am 2 t 4 u 4 - st 3 u 3 (7s 2 - Mu)) + 8m 8 s 7 - 192m 4 t 4 u 4 (2m 2 - s) - 8st 4 u 4 (s 2 + 3tu))/(D 2 TUt 2 u 2 ), 

f c o c o = (3 2 s(D2(16m 6 s 2 + 2m 4 s(3s 2 + 8tu) - m 2 s 2 (35s 2 - I06tu) + 8m 2 t 2 u 2 + s{t 2 + u 2 )(14s 2 - A7tu)) 
-3s(t-uf)/D 3 , 

f c o D o = -2f3 2 s(D 3 l6m 4 - D 2 2(Am 4 t(7t - 5u) - 2m 2 (7t 3 + U 2 u + 2tu 2 + u 3 ) + t(5t 3 - 23t 2 w + Atu 2 - 2u 3 )) 
-D2t(A8m 4 t 2 u + A3m 2 t 4 + m 2 u(7t 3 - \2t 2 u + Yitu 2 - 3u 3 ) + 14t 5 ) + 3t 3 (t - uf)/D 3 , 



f D -i = 8(2m 2 - s)(~DTU(6Am 6 (t 2 + u 2 ) - 8m 4 (5t 3 + 3t 2 u - 15tu 2 - u 3 ) + Am 2 t(3t 3 + 8tu 2 - 3u 3 ) 
+t 2 u(2t 2 - 5tu-8u 2 )) 

+TUu 2 (32m 2 t 2 {Am 4 + tz t ) - m 2 t 2 u{13t - 3u) + 2t 4 {t + Au)) - t 2 u 2 {Tu 3 {m 2 - 2t) + 2t 4 U))/(DTUt 2 u 2 ) 1 
f D -i B „ = 8(2m 2 - s)(D2m 2 (2m 4 {M - u) + \\m 2 t 2 - 2st 2 ) + TAm 2 t 2 (2st - u 2 ) + st 3 {3Tu + t 2 ))/{DTt 2 ), 



f D -i B o = 8(2m 2 - s)(-D(Am b (t + 5u) + 2m 4 (3r + 15iu + u 2 ) + m 2 t{3t 2 + u 2 ) - t 2 u{3t + 2u)) 

1m 1 

-t 4 (4m 4 -tu))/(DTt 2 ), 



f D -i Bl = 8(2m 2 - s)(DT2m 4 (8m 2 + 17t) - T2m 2 (m 4 u(t + 2u) + 7m 2 t 2 T + t 2 (t 2 - Au 2 )) + D3st 3 - Am 2 t 4 u)/(DTt 2 ), 



/n-i B i = 8(2m 2 - s)(Dm 2 (2m 2 (3t 2 + Atu + 2u 2 ) + t(At 2 + 3tu + 5u 2 )) 

lu 1 

+T(16m 8 s + 20m 6 st - 12m 4 t 2 u - Am 2 t 2 u 2 - t 3 u 2 ) + m 2 t 4 u)/(DTt 2 ) 



f D -i B o = -8(2m 2 - s){DA{m 2 s(i 2 + 2(m 2 s - t 2 )) + t(t - u) 3 )/(Df3 2 s), 
f D -i Bl =M2m 2 -s)z t /p 2 , 

f D -i c o = 8(2m 2 - s){D 2 Am 4 (3t - u) - D2t(m 4 (9t 2 - 2tu + u 2 ) - m 2 t{U 2 + lhtu + u 2 ) - t 2 {U 2 + 9tu + 2m 2 )) 
+8m 2 t 3 z t {2m 2 u - t 2 ) - t 5 {s 2 - Au 2 ))/(D 2 t), 

f D -i c o = -8(2m 2 - s)(D 2 Am 4 (t + 5u) - D2(m 4 u(9t 2 - 2tu + u 2 ) + m 2 t(2t 3 - 10t 2 u - Wtu 2 + u 3 ) + t 3 u(As - uj) 

lu 2 

+t 2 u(32m 4 uT - Am 2 t(2t - u)(t - u) - s 2 t 2 ))/(D 2 t), 
f D -^d = !6(2m 2 - s)(D(2m 6 (t 2 - Atu - u 2 ) + Am 4 t 2 (t - u) + m 2 t 2 {t - u) 2 - t 3 u 2 ) + f3 2 st 4 (Am 4 - tu))/(D 2 t), 
f D -i cl = 16(2m 2 - s)(D(2m 6 (3t 2 + u 2 ) + 2m 4 (t 3 - u 3 ) + m 2 t(2s(t 2 + u 2 ) - tu 2 ) - t 4 u) 

lu 2 

+z t t 3 u(m 4 + u 2 ) - 3m 4 st 3 u(3 2 + ra 2 tuT(2m 2 u 2 + t 3 ))/(D 2 t), 



f D -i c a = -A(2m 2 - s)(DA(2m 4 t(3t - 5u) + m 2 (3t 3 - 2t 2 u + 3tu 2 + 2u 3 ) + t(5t 3 + 2t 2 u + 2tu 2 - 3u 3 )) 
+8m 2 t 2 u((2m 2 - s)(t - 3u) + 2u 2 ) - t{t - uf)/D 2 , 

f D -i c i = 8(2m 2 - s){D(8m 6 s - 8m 4 t(3t - 2u) - m 2 (At 3 - 13t 2 u - \2tu 2 - 5m 3 ) - stu 2 ) 

1 5 

+2mV(10t 2 + u 2 ) - 10m 2 t 3 uT + s 3 t 2 u)/D 2 , 

f D -i c ° = 4 ( 2m2 - s){D(8m 4 (6t 2 - 5tu + u 2 ) - 2m 2 t(29t 2 - 15tu - 3u 2 ) - At(5t 3 + t 2 u - tu 2 - bu 3 )) 
+2m 2 t(t 3 T - 31u 3 T - 5t 3 u + 15tV + 6tu 3 - u 4 ) + t(t - u) 5 )/D 2 , 

f D -^d = -8(2m 2 - s)l3 2 s{D 2 2m 2 - D(m 2 (su + At 2 - Stu) - t 2 u) + (3 2 st 3 u - m 2 tu(t - u) 2 )/D 2 , 

f D -i D o = At(2m 2 - s)(-D(Am 4 (7t 2 - 2u 2 ) - 2m 2 (28t 3 + 23t 2 u + Wtu 2 + 2u 3 ) - 2t(6t 3 + \\t 2 u + 8tu 2 + 2w 3 )) 
+Am 4 tuT(9t - 7u) - 10m 2 t 4 T + lA m 2 t 2 u 2 (m 2 - t) + t 3 (s 3 + 2u 2 (s - t)))/D 2 , 

f D -i D o = At(2m 2 - s)(D(l2m 4 u(3t - u) + 2m 2 (2t 3 - 27t 2 u - Atu 2 - u 3 ) - 18t 3 u) - tu 2 z t (32m 4 - 17tu) 
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-4st 2 u 3 f3 2 - 2tu 2 T{3t 2 - 4u 2 ) + t 2 u(t 3 - t 2 u - u 3 ))/D 2 , 

f D -i D i = -8(2m 2 - s)(D 2 l6m 6 - Dm 2 (4m 4 (llt 2 + 5tu + 2u 2 ) + 4m 2 t(3t 2 + tu + 2u 2 ) - i(4i 3 + 3tu 2 + u 3 )) 
+mH 2 z t (9t 2 + 7u 2 ) - Am 2 t 3 u 2 T + st 2 u(s 2 T + tu(m 2 - t)))/D 2 , 

f D -i D x = -8(2m 2 - s)(D 2 l6m 6 - D(4m 6 (6t 2 + tu + 3u 2 ) + 2mV(5t + u) - m 2 t(2t 3 + 7t 2 u - u 3 )) 
-4m 4 t 3 uT + Um 4 t 2 u 2 (2m 2 - s) - Am 2 t 2 u 3 z t - 2m 2 t 5 u - s 2 t 3 u 2 )/D 2 1 

f D o = 4(D 2 T(512m 12 (t 2 + u 2 ) + 256m 10 i(4i 2 + 3tu + 7u 2 ) + 32m 8 (8i 4 + 37t 3 u + 67iV + 9tu 3 + 21u 4 )) 

-D 2 (8m 8 (7t 5 + 126t 4 u - 126iV - 2At 2 u 3 - 15tu 4 - 22u 5 ) + 4m 6 i(2i 5 - 2t 4 u - 8t 3 u 2 - 5t 2 u 3 - 63tu 4 - 74u 5 ) 
+m 4 2tu(m 5 + 433i 4 w - 291iV - 19t 2 u 3 - 6tu 4 - 4u 5 ) + m 2 2t 3 u(7t 4 + 70t 3 u + 52tu 3 + 2u 4 ) 
+t 3 u 2 (10t 4 + 30t 3 u + 37iV + 24tu 3 + 7w 4 )) 

+D0 2 sAm & t{128m 4 u 4 + m 2 tu(9lt 3 + 243u 3 ) - 2t 2 (2t 4 - 2t 3 u - 49u 4 )) 

-Dm 2 u 2 {lQ2Am 10 Tu 2 + 8m 4 tu(5t 4 + 3u 4 ) - 2m 2 t 5 (75t 2 - 61tu + M6uU + 522u 2 ) 
-t 3 (9t 5 + 75t 3 u 2 - lliV + 64tu 4 + 8u 5 )) 

-(3 2 m 2 st 3 u 4 (l96m 6 u - 472m 4 i 2 - 29m 2 i 2 u - 24i 4 - 33iV) 

+m 2 u 4 {16m 6 Tu 4 - lUm 2 t 7 - t 5 (184t 2 U - Atu 2 + 64Tu 2 - 49tt 3 )))/(Z) 2 ^ 2 sT[/tV), 

fD°D° = -(D 2 (128m w s + 128m 8 (3si - u 2 ) + Am 4 t{ht 3 + Aht 2 u - 63tu 2 - 3u 3 ) + 4m 2 t 3 (28t 2 + Altu + 26u 2 ) 
+2t 2 (l9t 4 + 2At 3 u + 15*V + 6tu 3 + u 4 )) 
+D(16m s (t(t 2 + u 2 )(29t + Uu) + 2u 2 (s 2 - tu)) - 52m 6 t(t 4 - u 4 ) - 4m 4 t 3 (66t 3 + 41i 2 u + tu 2 + 66u 3 ) 

-4m 2 i 3 (23< 4 + 5m 4 ) - 16< 8 ) + 66m 4 t e (t 2 + u 2 ) + 60m 2 t e Tu 2 + 24m 2 t e (t 3 + u 3 ) + 3t 6 {t 4 + u 4 ))/D 3 , 

fD°D° = -2(D 3 128m 6 {m 2 - 2s) + D 2 2m 2 (128m e tu + 8m 4 s(6t 2 + 6u 2 - tu) - m 2 (9s 4 - 17i 4 - 17u 4 + 106i 2 u 2 ) 

+2(2t 5 + t 4 u + 27iV + 27t 2 w 3 + tu 4 + 2u 5 )) 
+D2m 2 (256m 6 t 2 u 2 + 48m 2 t 3 u 3 - t 2 u 2 (26t 3 + 9stu + 26u 3 )) 
-3t 2 u 2 (30m 4 tu(t 2 + u 2 ) - 2m 2 {t 4 T + u 4 U) - s 4 tu))/D 3 . 
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